############################################
# NOTE: The replication is conducted in the order that tables and figures are placed in paper.
############################################

if (!require("pacman")) {
  install.packages("pacman")
  library(pacman)
}

# Load required packages
pacman::p_load(
  RColorBrewer, aod, arsenal, bookdown, broom, corrplot, devtools, dplyr, estimatr,
  evtree, fBasics, fabricatr, gganimate, ggplot2, ggstatsplot, gifski, grf, haven,
  hrbrthemes, rdrobust, kableExtra, knitr, lmtest, modelsummary, psych, randomizr,
  readr, readstata13, readxl, remotes, sandwich, sjPlot, sjlabelled, sjmisc, stargazer,
  texreg, tibble, tidyr, wesanderson, LexisNexisTools, quanteda, quanteda.textplots,
  readtext, quanteda.textstats, pwr
)

############################################ ################################################
######################################## MAIN TEXT ##########################################
############################################ ################################################

# Set working directory
setwd("~/Dropbox/Replication for Tureci-Sahin/data")

# Data
df <- readRDS(file = "~/Dropbox/Replication for Tureci-Sahin/data/df.rds")

# Function
normalize <- function(x) {
  data_range <- max(x, na.rm = TRUE) - min(x, na.rm = TRUE)
  return((x - min(x, na.rm = TRUE)) / data_range)
}

############################################ ################################################
######################################## MAIN TEXT ##########################################
############################################ ################################################

############################################ ################################################
#### Figure 1: Google search trends. ####
############################################ ################################################

google <- read.csv("multiTimeline.csv", sep = ';')
google$Date <- as.Date(google$Date, "%d.%m.%Y")
ggplot(google,aes(x=Date,y=Value,colour=Search,group=Search)) + geom_line() + geom_vline(xintercept=as.numeric(as.Date("0020-11-09")), color =  "grey",
                 linetype = "dashed") + theme_bw()  +
      labs(title = "Google Search for Pfizer–BioNTech developers",
           subtitle = "Grey line indicates the announcement of Pfizer–BioNTech vaccine success on November 9",
           y = "Participants",
           x = "Day") + 
  theme(axis.text.x = element_text(angle = 90))
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/googletrend-1.png")

############################################ ################################################
#### Figure 2: Word clouds with word frequency above 20 based on newspaper content. ####
############################################ ################################################

# List DOCX files in the directory
my_files <- list.files(pattern = ".DOCX", path = getwd(), full.names = TRUE, ignore.case = TRUE)

# Rename and read LexisNexis files
report <- lnt_rename(x = my_files, report = TRUE)
LNToutput <- lnt_read(
  x = getwd(),
  convert_date = FALSE,
  start_keyword = "auto",
  end_keyword = "auto",
  length_keyword = "auto",
  author_keyword = "auto",
  exclude_lines = "^LOAD-DATE: |^UPDATE: |^GRAFIK: |^GRAPHIC: ",
  recursive = FALSE,
  remove_cover = TRUE
)

# Extract metadata and articles
meta_df <- LNToutput@meta
articles_df <- LNToutput@articles

# Remove specific articles
articles_df <- articles_df[-c(80, 81), ]

# Create a corpus and tokenize the text
corp_inaug <- corpus(articles_df, text_field = "Article")
toks_immig <- quanteda::tokens(corp_inaug)

# Select tokens and remove unwanted words
toks <- quanteda::tokens_select(
  toks_immig,
  c("doct*", "uğur", "ugur", "şahin", "sahin", "özlem", "ozlem", "tureci", "company*", "corona", "virus",
    "pfizer", "vaccine*", "türeci", "euronews", "version", "biontech", "first", "two", "cancer", "can", "percent*",
    "end", "new", "also", "now", "body", "end", "copyright", "rights", "length", "november", "gmt", "end", "document",
    "reserved", "words", "wednesday", "us", "many", "one", "time", "year*", "date", "vacc*", "pandemic*", "doses",
    "pharmaceutical", "doses", "said", "million*", "world*", "according", "moderna*", "load*", "approval", "say*",
    "mrna", "euro*", "university", "already", "found*", "covid*", "around", "even", "still", "like", "however",
    "possibl*", "come*", "came", "section", "compan*", "medicine*", "billion*", "federal*", "cell*", "medic*",
    "develop*", "just*", "since*", "curevac", "test", "work", "people", "strüngmann", "reser*", "group*", "day",
    "thomas", "board", "monday", "rna", "week*", "month*", "berlin", "study*", "based*", "usa", "protect*", "health*",
    "long", "start*", "human*", "messenger", "money*", "news*", "human*", "data*", "money*", "gmbh", "byline",
    "immune", "risk", "tumor", "individual*", "fight", "drug", "immune*", "ahin", "share", "biotech", "early", "fact",
    "hospital", "candidate", "conclud*", "helmut", "kfw", "expect*", "stein*", "final", "köhler", "dpa", "bio*",
    "september", "infect*", "so-called", "explain*", "call*", "thing*", "today*", "hexal", "corona*", "high*", "day*",
    "bnt*", "homburg", "france", "+", "toilet", "influenza", "get", "based", "later", "p", "almost", "well", "without",
    "four", "want", "last", "without", "take", "much", "able", "co", "next", "make", "based*", "based ", "know", 
    "always", "received", "far", "actually", "three", "become", "together", "pdf*", "co*"),
  selection = "remove", padding = FALSE
)

# Create a document-feature matrix (dfm) and remove stopwords
dfmat_gr <- quanteda::tokens(toks, remove_numbers = TRUE,
                             remove_punct = TRUE,
                             remove_separators = TRUE,
                             remove_symbols = TRUE,
                             split_hyphens = TRUE,
                             remove_url = TRUE) %>% 
  tokens_remove(stopwords("en")) %>%
  dfm()

# Replace and standardize specific terms
dfmat <- dfm_replace(dfmat_gr, pattern = c("german", "germany", "germans"), replacement = c("german", "german", "german"))
dfmat <- dfm_replace(dfmat, pattern = c("turkey", "turkish", "turk", "turks"), replacement = c("turkish", "turkish", "turkish", "turkish"))
dfmat <- dfm_replace(dfmat, pattern = c("migrants", "immigrants", "migrant", "immigrant", "family", "guest", "worker"), replacement = c("migrant", "migrant", "migrant", "migrant", "migrant", "migrant", "migrant"))
dfmat <- dfm_replace(dfmat, pattern = c("success", "successful"), replacement = c("success", "success"))

# Create and save the word cloud
set.seed(100)
png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/media-3.png")
textplot_wordcloud(dfmat, min_count = 20, random_order = FALSE, color = RColorBrewer::brewer.pal(4, "Paired"))
dev.off()

############################################ ################################################
#### Figure 3: Data collection around the Pfizer-BioNTech vaccine success announcement.####
############################################ ################################################

p <- ggplot(df, aes(x=dDatum)) + 
      geom_bar() +
      geom_vline(xintercept = as.numeric(as.Date("2020-11-09")),
                 color =  "grey",
                 linetype = "dashed") +
      labs(title = "Interviews per day before and after treatment",
           subtitle = "Grey line indicates the announcement of Pfizer–BioNTech vaccine success on November 9",
           y = "Participants",
           x = "Day") + theme_bw()
p
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/descriptives-1.png", plot = p)

############################################ ################################################
#### Figure 4: Pfizer-BioNTech vaccine success announcement and attitudes towards immigration.
############################################ ################################################

# Model 1: Immigration Restriction Attitudes
model1<-lm_robust(Immigration_Restrict~Treatment, data=df, `se_type` = "stata")
model2<-lm_robust(Immigration_Restrict~Treatment, fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")

model3<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")

model4<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster, data=df, `se_type` = "stata")
model5<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20+occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
models1 <- list(
  "Baseline" = model1,
  "Baseline + FE" = model2,
  "Baseline + FE + controls" = model3,
  "Baseline + FE + controls + Day-cl SE" = model4,
  "Baseline + FE + controls + Day-cl SE + weights" = model5
)

plot1 <- modelplot(models1, coef_omit = "day|educ|age|gender|year|german|day|(Intercept)|occupation") +
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Immigration opportunities should be made easier') +
  guides(color = guide_legend(reverse = TRUE)) +
  geom_vline(xintercept = 0, color = 'grey') +
  coord_flip()
plot1 # 
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/immigration-1.png", plot = plot1)

# Model 2: Crime and Immigration Attitudes 
model6<-lm_robust(Immigration~Treatment, data=df, `se_type` = "stata")
model7<-lm_robust(Immigration~Treatment, fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model8<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model9<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster, data=df, `se_type` = "stata")
model10<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

models2 <- list(
  "Baseline" = model6,
  "Baseline + FE" = model7,
  "Baseline + FE + controls" = model8,
  "Baseline + FE + controls + Day-cl SE" = model9,
  "Baseline + FE + controls + Day-cl SE + weights" = model10
)

plot2 <- modelplot(models2, coef_omit = "day|educ|age|gender|year|german|day|(Intercept)|occupation") +
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Crime is increasing due to immigrants') +
  guides(color = guide_legend(reverse = TRUE)) +
  geom_vline(xintercept = 0, color = 'grey') +
  coord_flip()
plot2
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/immigration-2.png", plot = plot2)

############################################ ################################################
#### Figure 5: Pfizer-BioNTech vaccine success announcement and perceptions regarding immigration and immigrants.
############################################ ################################################

df$Factual_Info_Immigration <- df$CE50458/100
df$Competence_Immigration = (df$CE50457-1)/6

model36<-lm_robust(Factual_Info_Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model37<-lm_robust(Competence_Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
models3 <- list(
  "Perceived crime in relation to immigrants" = model36,
  "Perceived self-knowledge on immigration" = model37)

plot3 <- modelplot(models3, coef_omit = "day|educ|age|gender|year|german|day|(Intercept)|occupation")+
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Change in perceptions')+
  guides(color = guide_legend(reverse = TRUE))  + geom_vline(xintercept = 0, color = 'grey')+ coord_flip()
plot3
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/immigration perceptions-1.png", plot = plot3)

############################################ ################################################
######################################## APPENDIX ###########################################
############################################ ################################################

############################################ ################################################
#### Figure 6: The distribution of newspaper articles by date. ####
############################################ ################################################

dat_inaug <- read.csv("Media content analysis-1.de.en.csv", sep = ';')
corp_inaug <- corpus(dat_inaug, text_field = "English")
dat_inaug$Date <- as.Date(dat_inaug$Date, format = "%d.%m.%y")

docid <- paste(dat_inaug$Leaning,
               dat_inaug$Date,
               sep = " ")
docnames(corp_inaug) <- docid
docvars(corp_inaug, field = "Date")
tokeninfo <- summary(corp_inaug)
tokeninfo$Date <- as.Date(tokeninfo$Date, format = "%d.%m.%y")
dat <- aggregate(Tokens ~ Date, data = tokeninfo, FUN = "mean")

if (require(ggplot2)) {
  p <- ggplot(data = dat, aes(x = Date, y = Tokens, group = 1)) +
    geom_line() +
    theme_bw() +
    scale_x_date(date_labels = "%d.%m.%y", date_minor_breaks = "1 day")
  
  # Save the plot as 'media-1.png'
  ggsave("~/Dropbox/Replication for Tureci-Sahin/out/media-1.png", plot = p)
}

############################################ ################################################
#### Figure 7: Frequency of positive and negative words across newspaper articles. ####
############################################ ################################################

data_dictionary_LSD2015_pos_neg <- data_dictionary_LSD2015[1:2]
toks_gov_lsd <- tokens_lookup(toks, dictionary = data_dictionary_LSD2015_pos_neg)
# create a document document-feature matrix and group it by day
dfmat_gov_lsd <- dfm(toks_gov_lsd)

png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/media-2.png")
matplot(dfmat_gov_lsd, type = "l", lty = 1, col = 1:2,
        ylab = "Frequency", xlab = "")
grid()
legend("topleft", col = 1:2, legend = colnames(dfmat_gov_lsd), lty = 1, bg = "white")
dev.off()

############################################ ################################################
#### Figure 8: Wordfish fitted scaling model of all newspaper articles. ####
############################################ ################################################

tmod_wf <- textmodel_wordfish(dfmat_gr, dir = c(6, 5))
png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/media-5.png", width = 2000, height = 1500, res = 300)
textplot_scale1d(tmod_wf)
dev.off()

############################################ ################################################
#### Figure 9: Newspaper articles difference matrix based on Euclidean distance. ####
############################################ ################################################

toks_inaug <- quanteda::tokens(corp_inaug)
dfmat_gr <- quanteda::tokens(toks_inaug, remove_numbers = TRUE,
                             remove_punct = TRUE,
                             remove_separators = TRUE,
                             remove_symbols = TRUE,
                             split_hyphens = FALSE,
                             remove_url = TRUE) %>%
  tokens_remove(stopwords("en")) %>%
  dfm()

tstat_dist <- as.dist(textstat_dist(dfmat_gr, margin = "documents"))
hc <- hclust(tstat_dist, method = "complete")

png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/media-4.png", width = 1000, height = 2000, res = 300)
plot(hc, hang = -1, cex = 0.65, main = "Clustering Dendrogram")
dev.off()

############################################ ################################################
#### Table 1: Balance test of covariates by treatment. ####
############################################ ################################################

df$gender <- as.factor(df$gender_20)
df$educ_school <- as.factor(df$educ_school_20)
df$educ_job <- as.factor(df$educ_job_20)
df$occupation <- as.factor(df$occupation_20)
df$german_citizenship <- as.factor(df$german_citizenship_20)
df$marital_status <- as.factor(df$marital_status_20)

df_desc = df %>% select(gender, age, educ_job, educ_school, occupation, marital_status, german_citizenship, Treatment)
summary(tableby(Treatment ~ ., data = df_desc), title = "Balance test of covariates by treatment", text="latex")

wave43 = read.dta13("ZA7592_v2-2-0.dta")
wave43[wave43 %in% c(-99, -98, -90, -80)] <- NA
wave43 = wave43 %>% dplyr::select(id_g, AA43043)
df_secondary <- merge(df,wave43,by="id_g")
df_secondary <- df_secondary %>%
  mutate(Party_Affiliation = case_when(
    AA43043 == 2 ~ "SPD",
    AA43043 == 3 ~ "FDP",
    AA43043 == 4 ~ "Die Grunen",
    AA43043 == 5 ~ "Die Linke",
    AA43043 == 6 ~ "Other",
    AA43043 == 7 ~ "Other",
    AA43043 == 8 ~ "Other",
    TRUE ~ "CDU/CSU"
  ))

df_desc2 = df_secondary %>% select(Party_Affiliation, Treatment)
summary(tableby(Treatment ~ ., data = df_desc2), title = "Balance test of covariates by treatment", text="latex")

############################################ ################################################
#### Table 2: Balance test of covariates by treatment. ####
############################################ ################################################

# NOTE: The above code produces the content for both Table 1 and Table 2. 
# Due to space limitations, they are divided into two separate tables in the paper.

############################################ ################################################
#### Table 3: Summary statistics of the dependent variables. ####
############################################ ################################################

df_desc = df %>% select(CE50456, CE50461, CE50468, CE50047, CE50458, CE50457, CE50100, CE50464)
stargazer(df_desc)

############################################ ################################################
#### Table 4: The announcement of the vaccine and attitudes towards immigration (Welcoming immigration). #### 
############################################ ################################################

texreg(list(model1, model2, model3, model4, model5), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = " The announcement of the vaccine and attitudes towards immigration (Welcoming immigration).", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation", label = "tab:estimates-immigration-restrict", custom.gof.rows = list("State FE" = c("No", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("No", "No", "Yes", "Yes", "Yes"), "Clustered SE" = c("No", "No", "No","Yes","Yes"),  "Weight" = c("No", "No", "No","No","Yes")))

############################################ ################################################
#### Table 5: The announcement of the vaccine and attitudes towards immigration (Immigrants increase crime). ####
############################################ ################################################

texreg(list(model6, model7, model8, model9, model10), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the Vaccine and Attitudes towards Immigration (Immigrants increase crime)", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation|german|german", label = "tab:estimates-immigration-crime", custom.gof.rows = list("State FE" = c("No", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("No", "No", "Yes", "Yes", "Yes"), "Clustered SE" = c("No", "No", "No","Yes","Yes"),  "Weight" = c("No", "No", "No","No","Yes")))

############################################ ################################################
#### Table 6: The announcement of the vaccine and attitudes towards immigration. ####
############################################ ################################################

df_exc <- subset(df_secondary, dDatum > as.Date("2020-11-01"))
df_exc2 <- subset(df_secondary, german_citizenship_20==1)

#exclude November 1st
modelpc1<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_exc, `se_type` = "stata")
modelpc2<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_exc, `se_type` = "stata")
#exclude non-germans
modelpc3<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_exc2, `se_type` = "stata")
modelpc4<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_exc2, `se_type` = "stata")

#assume non-treatment on day of announcement 
df_secondary$Treatment[df_secondary$dDatum<"2020-11-10"] = 0
df_secondary$Treatment[df_secondary$dDatum=="2020-11-10"] = 1
df_secondary$Treatment[df_secondary$dDatum>"2020-11-10"] = 1
modelpc5<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")
modelpc6<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                    +occupation_20+german_citizenship_20,
                    fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

texreg(list(modelpc1, modelpc2, modelpc3, modelpc4, modelpc5, modelpc6),
       custom.model.names = c("Welcome", "Crime", "Welcome", "Crime", "Welcome", "Crime"),
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the Vaccine and Attitudes towards Immigration (Welcoming immigration)",
       stars = c(0.001, 0.01, 0.05, 0.1),
       omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation",
       label = "tab:estimates-rc-exclude1",
       custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Covariates" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), 
                              "Weight" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

# back to the former treatment date
df_secondary$Treatment[df_secondary$dDatum<"2020-11-09"] = 0
df_secondary$Treatment[df_secondary$dDatum=="2020-11-09"] = 1
df_secondary$Treatment[df_secondary$dDatum>"2020-11-09"] = 1

############################################ ################################################
#### Table 7:The Announcement of the vaccine success and attitudes towards environment (Environment and climate protection should be more important). ####
############################################ ################################################

model17<-lm_robust(Environment_Restrict~Treatment2, data=df, `se_type` = "stata")
model18<-lm_robust(Environment_Restrict~Treatment2, fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model19<-lm_robust(Environment_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model20<-lm_robust(Environment_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster, data=df, `se_type` = "stata")
model21<-lm_robust(Environment_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model17, model18, model19, model20, model21), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the vaccine success and attitudes towards environment (Environment and climate protection should be more important).", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation|german|german", label = "tab:estimates-rc-environment", custom.gof.rows = list("State FE" = c("No", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("No", "No", "Yes", "Yes", "Yes"), "Clustered SE" = c("No", "No", "No","Yes","Yes"),  "Weight" = c("No", "No", "No","No","Yes")))

############################################ ################################################
#### Table 8: The Announcement of the vaccine success and attitudes towards finance (Finance should be more important). ####
############################################ ################################################

model22<-lm_robust(Finance_Restrict~Treatment2, data=df, `se_type` = "stata")
model23<-lm_robust(Finance_Restrict~Treatment2, fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model24<-lm_robust(Finance_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), data=df, `se_type` = "stata")
model25<-lm_robust(Finance_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster, data=df, `se_type` = "stata")
model26<-lm_robust(Finance_Restrict~Treatment2+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model22, model23, model24, model25, model26), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the vaccine success and attitudes towards finance (Finance should be more important).", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation|german|german", label = "tab:estimates-rc-finance", custom.gof.rows = list("State FE" = c("No", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("No", "No", "Yes", "Yes", "Yes"), "Clustered SE" = c("No", "No", "No","Yes","Yes"),  "Weight" = c("No", "No", "No","No","Yes")))

############################################ ################################################
#### Table 9: Announcement of the vaccine success and randomly selected placebo variables (Part 1). ####
############################################ ################################################

selected_vars <- c(
  "Treatment", "gender_20", "age", "educ_school_20", "educ_job_20", "occupation_20",
  "german_citizenship_20", "state_20", "cluster", "webal",
  "expCE47383", "rndCE47402", "rndCE47403", "rndCE47404", "rndCE47405", "rndCE47406",
  "expZF50205", "rndZF50205", "ZF50205", "ZF50206", "ZF50207", "ZF50208", "ZF50209",
  "rndZF50210", "ZF50210", "ZF50211", "ZF50212", "ZF50213", "ZF50214", "rndZF50215",
  "ZF50215", "ZF50216", "ZF50217", "ZF50218", "ZF50219", "rndZF50220", "ZF50220",
  "ZF50221", "ZF50222", "ZF50223", "ZF50224", "expZF50225", "ZF50225", "ZF50226",
  "ZF50227", "ZF50228", "ZF50229", "ZF50230", "ZF50231", "ZF50232", "ZF50233", "ZF50234",
  "ZF50235", "ZF50236", "ZF50237", "ZF50238", "ZF50239", "ZF50240", "expCE50452_1",
  "expCE50452_2", "CE50452", "CE50047", "CE50100", "CE50453", "CE50454", "CE50463",
  "CE50464", "CE50465", "CE50466", "CE50467", "CE50468", "CE50469", "CE50471", "CE50472",
  "CE50473", "CE50475", "ZJ50053", "ZJ50054", "ZJ50055", "QE50001", "QE50002", "QE50003",
  "QE50004", "QE50005", "QE50006", "QE50007"
)

filtered_df <- df_secondary[selected_vars]
numeric_vars <- sapply(filtered_df, is.numeric)
filtered_numeric_df <- filtered_df[, numeric_vars]
vars_to_exclude <- c("Treatment", "gender_20", "age", "educ_job_20", "educ_school_20", "occupation_20", "german_citizenship_20", "state_20", "cluster", "webal")
all_vars <- names(filtered_numeric_df)
vars_to_keep <- setdiff(all_vars, vars_to_exclude)
models <- list()
for (var in vars_to_keep) {
  formula <- paste(var, "~ Treatment + gender_20 + age + educ_job_20 + educ_school_20 + occupation_20 + german_citizenship_20")
    model <- lm_robust(as.formula(formula),
                     fixed_effects = ~ factor(state_20),
                     cluster = cluster,
                     weights = webal,
                     data = df_secondary,
                     se_type = "stata")
    models[[var]] <- model
}

set.seed(123)

num_variables_to_analyze <- 10
selected_variables <- sample(vars_to_keep, num_variables_to_analyze)
selected_models <- list()
for (var in selected_variables) {
  formula <- paste(var, "~ Treatment + gender_20 + age + educ_job_20 + educ_school_20 + occupation_20 + german_citizenship_20")
    model <- lm_robust(as.formula(formula),
                     fixed_effects = ~ factor(state_20),
                     cluster = cluster,
                     weights = webal,
                     data = df_secondary,
                     se_type = "stata")
  
  selected_models[[var]] <- model
}

texreg(selected_models, booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Announcement of the vaccine success and randomly selected placebo variables",
       stars = c(0.001, 0.01, 0.05),
       omit.coef = "age|day|gender|education|job|school|marital|number|occupation|german|internet|occupation",
       label = "tab:random-variables")

############################################ ################################################
#### Table 10: Announcement of the vaccine success and randomly selected placebo variables (Part 2). ####
############################################ ################################################

### NOTE: The output continues here due to space constraints in the previous Table.

############################################ ################################################
#### Table 11: Announcement of the vaccine success and perception of threat from COVID-19. ####
############################################ ################################################

df$Covid_Threat = (df$ZJ50054/10)
df$Covid_Likelihood = (df$ZJ50055/7)

model42<-lm_robust(Covid_Threat~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model43<-lm_robust(Covid_Likelihood~Treatment+gender_20+age+educ_job_20+educ_school_20 +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster, data=df, `se_type` = "stata")

texreg(list(model42, model43), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Table 11: Announcement of the vaccine success and perception of threat from COVID-19.", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|day|gender|education|job|school|marital|number|occupation|german|internet|occupation", label = "tab:mechanism-threat", custom.gof.rows = list("State FE" = c("Yes", "Yes"), "Covariates" = c("Yes", "Yes"), "Clustered SE" = c("Yes", "Yes"),  "Weight" = c("Yes", "Yes")))

############################################ ################################################
#### Table 12: Alternative bandwidths post-treatment and attitudes towards immigration. #### 
############################################ ################################################

df$Treatment7[df$dDatum<"2020-11-09"] = 0
df$Treatment7[df$dDatum=="2020-11-09"] = 0
df$Treatment7[df$dDatum>"2020-11-09" & df$dDatum<"2020-11-17"] = 1
df$Treatment7[df$dDatum=="2020-11-17"] = NA
df$Treatment7[df$dDatum>"2020-11-17"] = NA

df$Treatment14[df$dDatum<"2020-11-09"] = 0
df$Treatment14[df$dDatum=="2020-11-09"] = 0
df$Treatment14[df$dDatum>"2020-11-09" & df$dDatum<"2020-11-24"] = 1
df$Treatment14[df$dDatum=="2020-11-24"] = NA
df$Treatment14[df$dDatum>"2020-11-24"] = NA

model11<-lm_robust(Immigration_Restrict~Treatment7+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model12<-lm_robust(Immigration_Restrict~Treatment14+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model13<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model14<-lm_robust(Immigration~Treatment7+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model15<-lm_robust(Immigration~Treatment14+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model16<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
texreg(list(model11, model12, model3, model14, model15, model16), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Alternative bandwidths post-treatment and attitudes towards immigration.", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation", label = "tab:estimates-rc-bandwidths", custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),  "Weight" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

############################################ ################################################
#### Table 13: Altering treatment day and attitudes towards immigration. 
############################################ ################################################

df$Fake_Treatment1[df$dDatum<"2020-11-16"] = 0
df$Fake_Treatment1[df$dDatum=="2020-11-16"] = 0
df$Fake_Treatment1[df$dDatum>"2020-11-16"] = 1

df$Fake_Treatment2[df$dDatum<"2020-11-21"] = 0
df$Fake_Treatment2[df$dDatum=="2020-11-21"] = 0
df$Fake_Treatment2[df$dDatum>"2020-11-21"] = 1

df$Fake_Treatment3[df$dDatum<"2020-11-26"] = 0
df$Fake_Treatment3[df$dDatum=="2020-11-26"] = 0
df$Fake_Treatment3[df$dDatum>"2020-11-26"] = 1

model1<-lm_robust(Immigration_Restrict~Fake_Treatment1+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model2<-lm_robust(Immigration_Restrict~Fake_Treatment2+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model3<-lm_robust(Immigration_Restrict~Fake_Treatment3+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model4<-lm_robust(Immigration~Fake_Treatment1+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model5<-lm_robust(Immigration~Fake_Treatment2+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
model6<-lm_robust(Immigration~Fake_Treatment3+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model1, model2, model3, model4, model5, model6), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Altering treatment day and attitudes towards immigration.", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation", label = "tab:estimates-rc-timing", custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Covariates" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),  "Weight" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

############################################ ################################################
#### Table 14: Announcement of the vaccine success and attitudes towards immigration, HTE by party affiliation. ####
############################################ ################################################

model5<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20+Party_Affiliation,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")
model6<-lm_robust(Immigration_Restrict~Treatment*Party_Affiliation+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")
model7<-lm_robust(Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20+Party_Affiliation,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")
model8<-lm_robust(Immigration~Treatment*Party_Affiliation+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

texreg(list(model5, model6, model7, model8), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Announcement of the vaccine success and attitudes towards immigration, HTE by party affiliation.",
       stars = c(0.001, 0.01, 0.05, 0.1),
       omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation",
       label = "tab:estimates-immigration-restrict-party-affiliation")


############################################ ################################################
#### Table 15: Announcement of the vaccine success and attitudes towards immigration, HTE by West vs. East Germany. ####
############################################ ################################################

df_secondary$West <- ifelse(df_secondary$state_20 %in% c(1, 4, 5, 6, 7, 8, 9), 1, 0)
df_secondary$East <- ifelse(df_secondary$state_20 %in% c(11, 13, 14, 15, 16), 1, 0)

modelwest1<-lm_robust(Immigration_Restrict~Treatment*East+gender_20+age+educ_job_20+educ_school_20
                      +occupation_20+german_citizenship_20,
                      fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

modelwest2<-lm_robust(Immigration~Treatment*East+gender_20+age+educ_job_20+educ_school_20
                      +occupation_20+german_citizenship_20,
                      fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

texreg(list(modelwest1, modelwest2), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the Vaccine and Attitudes towards Immigration",
       stars = c(0.001, 0.01, 0.05, 0.1),
       omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation",
       label = "tab:estimates-immigration-restrict-east")

############################################ ################################################
#### Figure 10: RD plot for immigration variables. ####
############################################ ################################################

png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/rd1.png", width = 1000, height = 800, res = 150)
rdplot(y = df$Immigration_Restrict,
       x = df$day,
       c = 10,
       p = 1,
       ci = 95,
       binselect = "qsmv",
       kernel = "uniform",
       title = "RD Plot for Vaccine Success Announcement",
       y.label = "Immigration opportunities should be made easier",
       x.label = "November 1st to November 30th, 2020")
dev.off()

png(filename = "~/Dropbox/Replication for Tureci-Sahin/out/rd2.png", width = 1000, height = 800, res = 150)
rdplot(y = df$Immigration,
       x = df$day,
       c = 10,
       p = 1,
       ci = 95,
       binselect = "qsmv",
       kernel = "uniform",
       title = "RD Plot for Vaccine Success Announcement",
       y.label = "Crime is increasing due to immigrants",
       x.label = "November 1st to November 30th, 2020")
dev.off()

############################################ ################################################
#### Figure 11: Sensitivity analysis. ####
############################################ ################################################

n1 <- 1561  
n2 <- 443 
alpha <- 0.05
power <- 0.8

effect_size <- pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha, power = power)$d
effect_size

alpha_05 <- 0.05
alpha_10 <- 0.1
power_values <- seq(0.6, 0.9, by = 0.01)

effect_sizes_05 <- numeric(length(power_values))
effect_sizes_10 <- numeric(length(power_values))

for (i in seq_along(power_values)) {
  effect_sizes_05[i] <- pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha_05, power = power_values[i])$d
  effect_sizes_10[i] <- pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha_10, power = power_values[i])$d
}

data_for_plot <- data.frame(Power = rep(power_values, 2),
                            EffectSize = c(effect_sizes_05, effect_sizes_10),
                            Alpha = factor(rep(c("0.05", "0.1"), each = length(power_values))))

closest_05 <- which.min(abs(effect_sizes_05 - 0.15))
closest_10 <- which.min(abs(effect_sizes_10 - 0.15))

power = ggplot(data_for_plot, aes(x = Power, y = EffectSize, color = Alpha)) +
  geom_line() +
  geom_point(data = data_for_plot[data_for_plot$Power == power_values[closest_05] & data_for_plot$Alpha == "0.05",], color = "red", size = 2) +
  geom_point(data = data_for_plot[data_for_plot$Power == power_values[closest_10] & data_for_plot$Alpha == "0.1",], color = "blue", size = 2) +
  theme_minimal() +
  scale_x_continuous(limits = c(0.6, 0.9)) +
  scale_y_continuous(limits = c(NA, 0.18)) +
  labs(title = "Effect Size by Power",
       x = "Power",
       y = "Expected Effect Size (Cohen's d)",
       color = "alpha") +
  theme(legend.title = element_text(size = 8))
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/Effect size and power.jpeg", plot = power, width = 10, height = 8, dpi = 300, device = "jpeg")

############################################ ################################################
#### Table 16: Announcement of the vaccine success and attitudes towards immigration, HTE by media consumption. ####
############################################ ################################################

df_secondary$Internet_use="Weekly"
df_secondary$Internet_use[df_secondary$internet_usage_20==8]="Daily"
df_secondary$Internet_use[df_secondary$internet_usage_20==7]="Weekly"
df_secondary$Internet_use[df_secondary$internet_usage_20==6]="Weekly"

model6<-lm_robust(Immigration_Restrict~Treatment*Internet_use+gender_20+age+educ_job_20+educ_school_20+occupation_20+german_citizenship_20, cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")
model7<-lm_robust(Immigration~Treatment*Internet_use+gender_20+age+educ_job_20+educ_school_20+occupation_20+german_citizenship_20, cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

texreg(list(model6, model7), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "The Announcement of the Vaccine and Attitudes towards Immigration",
       stars = c(0.001, 0.01, 0.05, 0.1),
       omit.coef = "age|educ|day|gender|marital|number|occupation|german|occupation",
       label = "tab:estimates-immigration-restrict-internet")


############################################ ################################################
#### Table 17: Announcement of the vaccine success and perceptions on immigration. ####
############################################ ################################################

df$Factual_Info_Immigration <- df$CE50458/100
df$Competence_Immigration = (df$CE50457-1)/6

model36<-lm_robust(Factual_Info_Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model37<-lm_robust(Competence_Immigration~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model36, model37), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Announcement of the vaccine success and perceptions on immigration.", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|day|gender|education|job|school|marital|number|occupation|german|internet|occupation", label = "tab:mechanism-perception-imm", custom.gof.rows = list("State FE" = c("Yes", "Yes"), "Covariates" = c("Yes", "Yes"), "Clustered SE" = c("Yes", "Yes"),  "Weight" = c("Yes", "Yes")))

############################################ ################################################
#### Table 18: Announcement of the vaccine success and perceived self-knowledge on other policy areas. ####
############################################ ################################################

df$Competence_Finance = (df$CE50100-1)/6
df$Competence_Environment = (df$CE50464-1)/6

model38<-lm_robust(Competence_Finance~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model39<-lm_robust(Competence_Environment~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model38, model39), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Announcement of the vaccine success and perceived self-knowledge on other
policy areas.", stars = c(0.001, 0.01, 0.05), omit.coef = "age|day|gender|education|job|school|marital|number|occupation|german|internet|occupation", label = "tab:mechanism-perception-other", custom.gof.rows = list("State FE" = c("Yes", "Yes"), "Covariates" = c("Yes", "Yes"), "Clustered SE" = c("Yes", "Yes"),  "Weight" = c("Yes", "Yes")))

############################################ ################################################
#### Table 19: Announcement of the vaccine success and political satisfaction. #### 
############################################ ################################################

df$Satisfaction_CDU = normalize(as.numeric(df$CE50153))
df$Satisfaction_SPD = normalize(as.numeric(df$CE50154))
df$Satisfaction_Gruene = normalize(as.numeric(df$CE50155))
df$Satisfaction_Linke = normalize(as.numeric(df$CE50156))
df$Satisfaction_AfD = normalize(as.numeric(df$CE50312))
df$Satisfaction_FDP = normalize(as.numeric(df$CE50313))

modelst1<-lm_robust(Satisfaction_CDU~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelst2<-lm_robust(Satisfaction_SPD~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelst3<-lm_robust(Satisfaction_Gruene~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelst4<-lm_robust(Satisfaction_Linke~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelst5<-lm_robust(Satisfaction_AfD~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelst6<-lm_robust(Satisfaction_FDP~Treatment+gender_20+age+educ_job_20+educ_school_20
                   +occupation_20+german_citizenship_20,
                   fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(modelst1, modelst2, modelst3, modelst4, modelst5, modelst6), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Table 19: Announcement of the vaccine success and political satisfaction.", stars = c(0.001, 0.01, 0.05, 0.1), omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation",
       label = "tab:mechanism-political-satisfaction",
       custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", 'Yes', 'Yes'),
                              "Covariates" = c("Yes", "Yes", "Yes", "Yes", 'Yes', 'Yes'),
                              "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", 'Yes', 'Yes'),
                              "Weight" = c("Yes", "Yes", "Yes", "Yes", 'Yes', 'Yes')))

############################################ ################################################
#### Figure 12: Perceived importance of the three policy areas. ####
############################################ ################################################

modelimm<-lm_robust(Immigration_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelenv<-lm_robust(Environment_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
modelfinan<-lm_robust(Finance_Restrict~Treatment+gender_20+age+educ_job_20+educ_school_20
                  +occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

models_compare <- list(
  "Immigration" = modelimm,
  "Environment" = modelenv,
  "Finance" = modelfinan)

plot_compare <- modelplot(models_compare, coef_omit = "day|educ|age|gender|year|german|day|(Intercept)|occupation")+
  labs(x = 'Coefficients', 
       y = 'Term names',
       title = 'Perceived importance of policy areas')+
  guides(color = guide_legend(reverse = FALSE))  + geom_vline(xintercept = 0, color = 'grey')+ coord_flip()
plot_compare
ggsave("~/Dropbox/Replication for Tureci-Sahin/out/perceived importance.png", plot = plot_compare, width = 10, height = 8, dpi = 300)

############################################ ################################################
#### Table 20: Announcement of the vaccine success and PCA variable of attitudes towards immigration. #### 
############################################ ################################################
df_secondary$Factual_Info_Immigration <- df_secondary$CE50458/100
df_secondary$Competence_Immigration = (df_secondary$CE50457-1)/6
df_secondary$Flipped_Immigration <- 1 - as.numeric(df_secondary$Immigration)
df_secondary$Flipped_Competence_Immigration <- 1 - as.numeric(df_secondary$Competence_Immigration)
df_secondary$Flipped_Factual_Info_Immigration <- 1 - as.numeric(df_secondary$Factual_Info_Immigration)

immigration_vars <- c("Flipped_Immigration", "Immigration_Restrict", "Flipped_Competence_Immigration", 'Flipped_Factual_Info_Immigration')
complete_cases <- complete.cases(df_secondary[, immigration_vars])
pca_result <- prcomp(na.omit(df_secondary[immigration_vars]), scale. = TRUE)
pc1_scores <- pca_result$x[, 1]
pc1_scores_normal = normalize(pc1_scores)
summary(pc1_scores_normal)
df_secondary$PC1 <- NA
df_secondary[complete_cases, "PC1"] <- pc1_scores_normal*-1

model5<-lm_robust(PC1~Treatment+gender_20+age+educ_job_20+educ_school_20+occupation_20+german_citizenship_20,
                  fixed_effects = ~ factor(state_20), cluster=cluster,weights = webal, data=df_secondary, `se_type` = "stata")

texreg(list(model5), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Table 20: Announcement of the vaccine success and PCA variable of attitudes towards immigration.", stars = c(0.001, 0.01, 0.05, 0.1),
       omit.coef = "age|educ|day|gender|marital|number|occupation|german|internet|occupation|german|german",
       label = "tab:estimates-immigration-pca",
       custom.gof.rows = list("State FE" = c("Yes"),
                              "Covariates" = c("Yes"),
                              "Clustered SE" = c("Yes"),
                              "Weight" = c("Yes")))
